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Abstract 

Tissue microstructure, in particular the alignment of myocytes (fibre direction) and their lateral organisation into sheets, is 
fundamental to cardiac function. We studied the effect of microstructure on contraction in a computational model of rat left 
ventricular electromechanics. Different fibre models, globally rule-based or locally optimised to DT-MRI data, were 
compared, in order to understand whether a subject-specific fibre model would enhance the predictive power of our model 
with respect to the global ones. We also studied the impact of sheets on ventricular deformation by comparing: (a) a 
transversely isotropic versus an orthotropic material law and (b) a linear model with a bimodal model of sheet transmural 
variation. We estimated ejection fraction, wall thickening and base-to-apex shortening and compared them with measures 
from cine-MRI. We also evaluated Lagrangian strains as local metrics of cardiac deformation. Our results show that the 
subject-specific fibre model provides little improvement in the metric predictions with respect to global fibre models while 
material orthotropy allows closer agreement with measures than transverse isotropy. Nonetheless, the impact of sheets in 
our model is smaller than that of fibres. We conclude that further investigation of the modelling of sheet dynamics is 
necessary to fully understand the impact of tissue structure on cardiac deformation. 
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Introduction 

Cardiac tissue is arranged into a complex architecture. Due to 
their elongated shape, myocytes are highly anisotropic and their 
local orientation in cardiac tissue is usually called the fibre 
direction [1]. Furthermore, the extracellular connective tissue 
matrix that surrounds myocytes groups them into layers 4-6 cells 
thick [2], separated by cleavage planes, forming a branching 
laminar architecture. Such layers are commonly referred to as 
sheets. Although there is general agreement about the local 
alignment of cells and their organisation into sheets, disagreement 
remains on their exact characteristics in beating hearts, their 
alteration due to pathology and in general their degree of 
contribution to cardiac function [3]. Tissue microstructure 
contributes to the heterogeneity of tissue conductivity and 
therefore directly affects electrical wave propagation, with 
dramatic consequences in the case of pathology. Regarding 
mechanics, the importance of taking into account fibre and sheet 
orientation in the definition of passive and active behaviour of 
cardiac tissue has already been highlighted [4] . It is along the fibre 
direction that electrical excitation propagates fastest and that the 



principal component of contraction takes place. When included in 
an electromechanical model, fibres are often mathematically 
prescribed [5] . Recent advances in imaging technology have made 
it possible to measure local microstructure from data such as ex vivo 
diffusion tensor MRI (DT-MRI). The primary eigenvector of the 
diffusion tensor has been shown to correlate with the fibre 
direction in the tissue [6]. Electromechanical models based on 
DT-MRI data suggest that fibres are a dynamical system in 
cardiac structure, undergoing adaptation that minimises shear 
strains, behaviour that can have important consequences in tissue 
remodelling in pathology [7] . The importance of local myocyte 
orientation in cardiac function, both in health and disease, has 
promoted considerable effort in recent years towards subject- 
specific modelling of cardiac electromechanics [8], [9], [10]. 
Pathological conditions such as dyssynchrony [1 1] have been 
investigated, with realistic ventricular to whole-organ geometries, 
combined with fibre orientations obtained from imaging [10]. The 
results are promising, but it is an open question as to whether 
current electromechanical models are becoming a useful support 
tool in clinical decision-making. 
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The laminar structure of myocardial tissue is also important for 
cardiac electromechanics. It has been established that sheets 
undergo shortening in the fibre direction, extension in the in- 
sheet/cross-fibre plane, thickening and shear, whose regional 
variability and transmural distribution are important for the 
distribution of strain and stress throughout the ventricle during 
cardiac motion [12]. More recendy, changes in sheet intersection 
angles have been identified to occur during the cycle of 
contraction and relaxation [13]. These sheet dynamics are thought 
to form one of the underlying mechanisms to explain wall- 
thickening during contraction [14]. Sheets have been modelled as 
continuously varying transmurally [2], although more recent 
experimental findings suggest that the transmural distribution of 
sheets is more complex and discontinuous [13], see for example 
Figure 1(f). Sheet orientations directly obtained from histological 
data, usually from a small number of ventricular locations only, 
have been employed in computational studies of the effect of tissue 
microstructure on electrophysiology [15], but in general the 
importance of sheets has not been fully acknowledged in the 
cardiac modelling community in spite of the experimental 
evidence. 

Two motivations underlie our study. First, we aim to 
understand whether the use of subject-specific fibres would 
enhance the predictive power of our model, or whether it would 
be as informative as simple globally prescribed fibres within our 
modelling framework. To this purpose, globally defined rule-based 
fibres were compared with a subject-specific model of fibres that is 
locally optimised in order to minimise the discrepancy with DT- 
MRI data [5]. 

Our second motivation is the exploration of the often 
underestimated importance of sheet dynamics in the common 



electromechanical models. We want to test the impact of sheet 
orientation in the mechanics of the left ventricle. This was 
achieved by (a) comparing transverse isotropic with orthotropic 
material laws for cardiac tissue and (b) under the orthotropic 
material properties, testing two global models of transmural 
variation of the sheet orientation, namely a linear and a bimodal 
model. The aim of the latter comparison was to see whether the 
latest findings regarding the discontinuous multiple-population 
distribution of sheets [13] (mimicked by the bimodal model) can 
improve the accuracy of the electromechanical model, as opposed 
to the standard smooth variation of sheet orientation [2] (the linear 
model). 

In order to quantify the effect of the different models of fibre and 
sheet orientation on contraction, global and local metrics were 
considered. The global metrics are ejection fraction, wall thicken- 
ing, and shortening in the base-to-apex direction. We compared 
their predicted values with those obtained from rat cine-MRI data. 
The local metrics considered are the Lagrangian circumferential 
strain, as a measure of contraction, circumferential-radial shear 
strain, as a measure of the relative twist of sub-endocardial and sub- 
epicardial layers, and circumferential-longitudinal shear strain, as a 
measure of torsion [16]. 

A preliminary study of the effect of fibre direction and material 
orthotropy on contraction had been performed with a simplified 
electromechanical model (without pressure boundary conditions) 
[17]. The current study represents a significant extension in our 
model development, to further the study of the impact of sheets 
and the analysis of the local effect of tissue structure on 
deformation. The contribution of the study is twofold. First, we 
show that the fibre model affects significantly the prediction of the 
metrics considered. Therefore, it is important to choose the most 




16mm 




wdp=n 




Figure 1. Modelling ventricular geometry and tissue structure, (a) The ellipsoidal mesh seen in a transmural-longitudinal cut. (b) Basal slice of 
the mesh cut parallel to the short axis, (c) Midventricular slice from rat cine-MRI at beginning of contraction, (d)-(e) Helix, and sheet angle definition, 
(f) Image of transmural histological sample highlighting the sheets. 
doi:1 0.1 371 /journal.pone.0092792.g001 
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suitable fibre setting. On the other hand, our results question the 
efficacy of using personalised myocyte orientations when com- 
pared to the limitations of state-of-the-art electromechanical 
models, since they do not seem to provide substantial improve- 
ments, both in terms of global and local measures, compared to 
simpler global rule-based fibre definitions, at least for the healthy 
left ventricle of rat. Secondly, we highlight the need for further 
development of models of sheets. On the one hand, the 
introduction of sheet orientation in the model by switching from 
transverse isotropy to orthotropy shows a small but consistent 
improvement in the metrics predictions. On the other hand, we 
found that, although experimentally more realistic, the bimodal 
model for sheet orientation did not improve and in some cases 
counteracted the effect of introducing material orthotropy. We 
believe that such an important structural change should affect 
more evidently the predicted deformation patterns. Given the 
results presented in this study, we conclude that our current 
electromechanical modelling approach may not be able to 
successfully make use of subject-specific tissue structure, with all 
its additional complexity and computational cost, unless a better 
understanding of the biological roles, and the computational 
modelling of sheets and their mechanical properties is achieved. 

Materials and Methods 

No new animal experiment was required for the results 
presented in this paper. The images used come from a dataset 
previously published and for which the ethics statement was duly 
reported [18]. However, with respect to those experiments, the 
authors confirm that all animal work was conducted in accordance 
with Schedule 1 of the UK Home Office Guidance on the 
Operation of Animals (Scientific Procedures) Act of 1986, and was 
approved by the Oxford University ethical review board. 

MRI dataset 

In vivo cine-MRI and subsequent ex vivo DT-MRI of rat hearts 
were acquired on a 9.4T (400 MHz) MR system (Agilent, Santa 
Clara, CA, USA). The cine-MRI dataset used a 2D multi-frame 
gradient echo sequence, and consisted of a stack of eight 
contiguous short-axis slices covering the heart from apex to the 
atrioventricular separation. Each short-axis slice comprised of 32- 
36 frames acquired throughout the cardiac cycle with an in-plane 
resolution of 100x100 fim and a slice thickness of 1.5mm. A short- 
axis view in end-diastole is shown in Figure 1(c). The DT-MRI 
images where subsequently acquired on excised rat hearts, fixed 
either in resting state or in contracture and embedded in agarose, 
at an isotropic resolution of 100 ^tm as described in [18]. 

Geometrical mesh for rat left ventricle 

We approximate the left ventricle with a thick-walled truncated 
ellipsoid. A first coarse tetrahedral mesh was generated with semi- 
axes and wall thickness estimated from the end-diastolic frame of 
one of the cine-MRI rat scans described (Figure 1(c)), which we 
took as the reference mesh (Figures 1(a) -1(b)). In order to estimate 
the shape of a mesh representative of the left ventricle in absence 
of blood pressure (called pre-P mesh), an iterative approach was 
applied. The reference mesh was shrunk by applying a uniform 
scaling matrix to give the pre-P mesh, which was then allowed to 
expand passively under a constant pressure of — 1 KPa applied to 
the endocardial surface to give the post-P mesh (the negative sign 
means that the pressure is applied inwards along the normal to the 
surface). This step was repeated for several iterations adjusting the 
scaling factor until we obtained an optimal pre-P mesh producing 
a post-P mechanics mesh close to the reference mesh. Details 



about the pre-P and post-P meshes are given in Table 1. The 
mechanics mesh for the electromechanical simulations is the pre-P 
mesh resulting from the optimisation performed, and has 
quadratic elements. The mesh for the electrophysiological 
component is instead obtained by refining the post-P mechanics 
mesh and has approximately 1.6-10 4 linear tetrahedral elements 
with a mean edge length of 5J-\0~ 2 cm. All meshes were 
generated with the open-source 3D finite element grid generator 
Gmsh [19]. 

Tissue structure definition 

Tissue structure is defined at the centroid of each element of the 
mechanical mesh by assigning fibre, sheet and sheet-normal 
directions. This orthonormal local system of coordinates is 
assigned according to a well-established approach [2]. First, the 
orthonormal system of material coordinates (radial u, circumfer- 
ential v and longitudinal w directions) for the centroid of a given 
element is identified. A first rotation about the radial direction u is 
performed, that transforms the vectors v and w into v a and w a . 
The rotation angle a is the helix angle that defines the fibre direction 
f = V;,. The second step consists of a rotation about the fibre 
direction f of the angle ft, the sheet angle. The vector u is rotated to 
Ub that represents the sheet direction s, while the vector w a is 
rotated into w a) s that gives the sheet-normal direction n. 
Figures 1(d) and 1(e) illustrate the first and second step, 
respectively. 

Different models for the fibre and sheet orientation are 
introduced in this tissue structure generation framework by 
modelling the transmural variation of the helix and sheet angle 
as described below. 

Models of fibre direction. A set of rule-based globally 
defined fibres were generated by assigning to each element of the 
mesh a value of the helix angle cc, defined as the function 
oc = oc(R,n,d) and calculated by eqn. (1) 

ot = R*sgn(l-2d)\\-2d\" (1) 

where sgn stands for the sign function, d is the normalised 
transmural depth, varying from 0 to 1 from endocardium to 
epicardium, R is the transmural rotation, i.e. the maximum 
absolute value of a allowed by the model, which will span between 

Table 1. Mechanical mesh properties. 



pre-P mesh post-P mesh 



Chpt 


0.34 


0.40 


bepi 


0.48 


0.48 




0.95 


0.99 




0.16 


0.26 


bf„do 


0.27 


0.3 




0.79 


0.84 


mean wall thickness 


0.18 


0.16 


Cavity volume (cm 3 ) 


0.07 


0.132 


Tissue volume (cm 3 ) 


0.267 


0.267 


Mean edge length 


7.2T0~ 2 


7.1T0- 2 


No. of elements 


8225 


8225 



Mechanical mesh properties before (pre-P) and after (post-P) passive expansion 
under pressure. Length unit is cm. See Figures 1(a) -1(b) for reference to the 
ellipsoid parameters. 
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— R at epicardium and +R and endocardium. Finally, the 
exponent n modulates the transmural variation of a, so for n = 1 it 
will be linear, n = 3 cubic. This is a standard approach for globally 
prescribed fibres introduced by Streeter et al. [1]. Four values of R 
were tested, R = 50, 60, 70, 80 degrees, spanning the range of 
values found experimentally [5] . For each value of R, four values 
of the exponent n were considered, « = 0.5, 1,2,3, for a total of 16 
globally defined fibre models. In the rest of the paper, these models 
are referred to as SnxRy, where S stands for Streeter-based and 
nxRy refers to the value of n and R used. For example, the fibre 
model Sn0.5R70 is obtained by imposing R = 10 and n = 0.5 in 
eqn. (1). Figure 2(a) shows the transmural variation of a. for a 
subset of these fibre models. 

The global fibre models are compared with a subject-specific 
fibre model developed by Karadag et al. [5]. This is based on a 
local optimisation of R and n carried out on 6 subjects from the 
DT-MRI rat dataset previously described. The cine-MRI rat scan 
used to generate the mesh comes from one of the subjects from this 
set. The left ventricle was divided into 1 6 regions, an adaptation of 
the 17-region classification of the American Heart Association 
[20], and on each of them eqn. (1) was fitted in order to minimise 
the angular difference with the data. The result is a 16-regions 
model of fibres angles, with a single couple (R,n) obtained for each 
region. We refer to this model as Karadag. 

Models of sheet orientation. The standard approach to the 
modelling of sheet orientation is that of a smooth variation in the 
transmural direction [2]. Recent experiments question this 
modelling approach showing the presence in the tissue of sheet 
populations with opposing orientation angles [3] and a strong level 
of discontinuity [13]. In order to compare the effect of the two 
different approaches, we have compared a simple linear transmu- 
ral variation of the sheet angle /? with a bimodal model. The linear 
model for /? goes from B= —80 at endocardium to B= +80 at 
epicardium, and is described by P = B{2d—\). The bimodal 
model consists of a continuous approximation obtained by means 
of the sigmoidal function in eqn. (2) 

P = B{2H-\) where H = (\ +exp(-K(2d- l)))" 1 (2) 



where K is a parameter affecting the steepness of the sigmoidal 
function. We chose K = 20 in order to have a distinct but smooth 
step in the trace. Figure 2(b) shows the transmural variation of the 
sheet angle /? for the two models used in this study. In the 
following, unless it is specified, the linear model for the sheet angle 
is used. 

Electromechanical model 

Model of electrophysiology. At the tissue level the mono- 
domain equations (eqn. (3) and (4)) for electrical wave propagation 
are solved, 

Z(Qf +h,n(g,V))-S/iaWV) + I slim = 0 (3) 

§=f(g,^) (4) 

where V is the transmembrane voltage [m V), t is time (m s), C m is 
the membrane capacitance (fiF/cm 2 ), x is the surface area-to- 
volume ratio (cm -1 ), a is the effective conductivity (mS/cm), a 
matrix encoding the conductivity along the fibre, sheet and sheet- 
normal directions. The function Iion(&V) is the ionic current 
provided by the cell model (piA/cm 2 ) and I s a m is the external 
stimulus (p.Ajcrrv'). The vector g represents the set of state 
variables forming the cell electrophysiological model, whose 
temporal variation is regulated by the ordinary differential 
equation (ODE) system in eqn. (4). The model of cellular 
electrophysiology is that of Wang and Sobie [21]. Table 2 
provides the parameters values we used for eqn. (3). The value of 
the surface area-to-volume ratio x is typical for rat myocytes [22]. 
The conductivities were adjusted in order to ensure that the 
resulting conduction velocity in the fibre and cross-fibre directions 
was within experimental range, with ratio — — — «2.3 [23]. 

cv cross— fibre 

Cubic wedge simulations of electrophysiology with straight fibres 
parallel to one of the coordinates axes were used for this purpose. 



A B 




0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 

transmural d (endo=0, epi=1 ) transmural d (endo=0, epi=1 ) 



Figure 2. Helix and sheet angle variation, (a) Transmural models of the helix angle tz for four of the global rule43ased fibre models, (b) 
Transmural models of sheet angle /}. 
doi:1 0.1 371 /journal.pone.0092792.g002 
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Table 2. Parameters of the monodomain equation (eq. (3)). 




Parameter 


Value(s) 




l.OriF/cm 2 


X 




/v/7111 


6.0' lO 5 /^/™ 3 


stim duration 


I ms 


Of 


3.0mS/cm 




l.OmS/cm 


a„ 


l.OmS/cm 


doi:1 0.1 371 /joumal.pone.0092792.t002 



Model of mechanical behaviour. The mechanical compo- 
nent is coupled to the electrical component via the transmembrane 
voltage. At the tissue level, the mechanical behaviour is modelled 
according to finite elasticity theory. Let us consider two systems of 
Cartesian coordinates, a reference system defined on the 
undeformed body, X = (Xi,X2,Xi), and a second system repre- 
senting the deformed state, x = (xi,X2,x^). The components of the 
deformation gradient tensor F, which describes the transformation from 
the undeformed state to deformed one, can be written as 

-f/M = "5 — ~- The Green deformation tensor C = F T F, i.e. 

Cmn= ^2iFiAiFiN, describes the deformation in the undeformed 
coordinate system. Finally, the Lagrangian strain tensor is 

E= — (C— I), where / is the identity matrix. 

The 2 nd Piola-Kirchoff stress tensor T for the total stress 
produced by the tissue is the sum T = TP" ssh ' e + T" c,h ' e of the stress 
due to the tissue passive response, T passn ' e , and active contraction, 
T ac " re '. The constitutive relation between passive stress and strain 
dW 

is given by T pasm ' e = ~^r> where W = W(E) is the exponential 
strain-energy function used by Usyk et al. [4] and given in eqn. (5), 

W(E)= ( ^-(e Q -\) + C h {J\n{J)-J+\) where (5) 

Q= EE b «4 {iJ}=UM and h=h 

> J 

where J = det(F), C a = 0.88kPa is a stiffness constant and h» for 
{ij} = {f,s,n} are stiffness parameters. The parameter Cj 
regulates the level of compressibility allowed by the material law: 
the higher the value the closer it is to being incompressible. After a 
preliminary analysis of the effect Q, on the overall compressibility 
of the tissue we chose Ci, = 200 kPa, since it ensures a negligible 
volume reduction VR=[(,V md -V start )/V start }100%= - 1.63% 
between beginning and end of contraction. Furthermore the 
distribution of the Jacobian of the deformation / over the mesh 
elements between beginning and end of contraction has an 
average 5 th — 95 percentile range of 0.95-1 .00, confirming that 
the model is nearly incompressible. 

Two models of material passive behaviour were tested 
representing transversely isotropic and orthotropic behaviour. 
The stiffness parameters for the orthotropic case are taken from 
[4] . The parameters for the transversely isotropic law are obtained 



by averaging the relevant parameters of the orthotropic law. The 
parameter values for both laws are reported in Table 3. 

The active stress at the tissue level, T ac,m ', depends on the active 
tension generated at the cell level. The cell contraction model we 
used was originally developed by Kerkchoffs [24] for canine 
myocytes. The model parameters associated with the mechanical 
activation time in the twitch component of the model, [b,t r and tf), 
were reduced by 50% in order to get a shorter contraction 
duration (as 150 ms), more compatible with the rat cardiac cycle 
frequency, while still generating an analogous amount of active 
tension ( » 1 6kPa at peak contraction) in that shorter interval. 
Finally, the total stress T is required to satisfy the conservation of 
momentum as shown in eq. (6) 

E ^(T MN F lN ) = 0 1=1,2,3 (6) 

N,M= 1,2,3 0Am 

where the effect of inertia and gravity is neglected. The system is 
solved with appropriate Dirichlet boundary conditions for 
displacement and Neumann boundary conditions for tractions, 
in this case the blood cavity pressure onto the endocardial surface, 
which is defined as follows. 

Cavity pressure trace. The temporal evolution of the intra- 
cavital blood pressure is approximated by means of a sigmoidal 
function ranging from — 1 kPa to — 5 kPa, where the negative 
sign is due to the fact the pressure is applied in the inward 
direction along the normal to the surface. The range of pressure 
applied is significantly lower than that typically measured in the 
rat left ventricle, which reaches approximately 16 kPa (or more, 
for example during physical activity) at end-systole [25] . This is a 
known limitation of some contraction models and follows existing 
computational modelling literature [26]. However, the chosen 
pressure trace ensures that the left ventricle undergoes first an 
isovolumic contraction and then an ejection phase. Figure 3(a) 
shows the pressure trace with an example of the resulting cavity 
volume. 

Electromechanical simulations 

All simulations were run using Chaste, an open-source 
simulation package for computational biology developed in our 
group (http://www.cs.ox.ac.uk/chaste) [27]. The simulation 
duration is 150 ms which is approximately the time at which full 
contraction is achieved. The timesteps used for solving the 
electrophysiology are 5-10~ 3 /*ra for the cell model (eq. (4)) and 
1-10 ms for the monodomain equation, while for the mechanics 
they are 0. 1 ms for the contraction model and 1 ms for eq. (6). The 
numerical scheme used has been described in previous work [28]. 



Table 3. Stiffness parameters values for the material law (eq. 
(5)). 




Parameter 


transverse Isotropy 


orthotropy 


% 


6.0 


6.0 


bss 


4.3 


7.0 


b, m 


4.3 


3.0 


b fs 


7.5 


12.0 


bf„ 


7.5 


3.0 


b„ 


4.3 


3.0 



doi:1 0.1 371 /joumal.pone.0092792.t003 
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time (ms) 



B 

62 r 
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maximum absolute R angle 



Figure 3. Pressure and volume traces and global metric EF. (a) Pressure trace and resulting cavity volume trace for simulation with fibre 
Sn1R70 and transverse isotropy. (b) EF for the different fibre models, with transversely isotropic material law. 
doi:1 0.1 371 /journal.pone.0092792.g003 



A single electrical stimulus of 6-10 5 f.iAcm~ 3 intensity is applied to 
the apical region, identified as the lower 5% of the left ventricle. 
This type of activation region is commonly used in experiments 
with Langendorff-perfused hearts, and it was chosen to mimic the 
apex-to-base direction of the electrical activation in the left 
ventricle. Dirichlet boundary conditions fix the nodes at the base 
of the ventricle in the long-axis direction, so that they are allowed 
to move in the short-axis plane. 

Global and local metrics 

The global measures used to evaluate the simulations are: 
ejection fraction EF, wall thickening WT and shortening in the 
base-to-apex direction SL. All the quantities are computed by 
comparing the state of the mesh at the beginning (to = 0 ms) and 
end (//■ = 150 ms) of contraction. EF is calculated as the relative 
difference between the mesh cavity volume vol at to and tf, 
EF= [(vo/, 0 — vol, f )/vol, 0 ] 100%. WT is the difference in wall- 
thickness W, WT =[(w tf — w lt> )/wt 0 ] WOVo, where W at a given 
time is obtained by averaging four measures of wall-thickness 
taken symmetrically onto a short-axis slice in the midventricular 
region (at approximately 25% of the long axis from the base). We 
chose this specific short-axis slice to compare with the cine images, 
where WT is more easily measured in the midventricular region. 
SL is calculated by comparing the z coordinate of the epicardial 
apical node along the z axis (which coincides with the long axis of 
the left ventricle, see Figure 1(a)) at to and tf, 
SL= [(z lf — Z( 0 )/z^] 100%. This is an adequate measure of the 
shortening along the base-to-apex direction since the base is kept 
fixed in the model and the apex moves almost exclusively along the 
long axis (of course, in situ, the movement is in the opposite 
direction, i.e. the atrio-ventricular valve plane moves 'downwards' 
in the direction of the apex, but in terms of changes in apico-basal 
length, this parameter offers a satisfactory approximation within 
the modelling framework used). Here SL<0 represents shorten- 
ing, while SL>0 represents elongation of the ventricle during 
contraction. 

In order to compare the predicted values with experimental 
data, we measured these metrics on the cine-MRI dataset used for 
building the ellipsoidal mesh by applying analogous methods. The 



measured values are EF cine = ll%, WT cine = 5A% which are 
consistent with values in the literature [29], and SL c j ne = — 10%. 

As local metrics we considered the Lagrangian circumferential 
strain Ecc, and two shear strains, the circumferential-longitudinal 
shear strain Eel, and the circumferential-radial shear strain Ecr. 
The Lagrangian strain tensor E in the Cartesian system is 
obtained from the deformation tensor F calculated during the 
simulation, that was sampled every 10ms. Then the conversion to 
cylindrical coordinates was applied. We could not compare these 
strain predictions directly with data, but they show good 
agreement with results in the literature as further explained in 
the Discussion. 

Results 

In this section we present the results of the series of 
electromechanical simulations we performed in order to (a) 
compare sixteen globally defined fibre models (SnxRy for 
x = 0.5,l,2,3 and y = 50,60,70,80) against one subject-specific 
fibre model (Karadag), and (b) evaluate the effect of introducing 
sheets by comparing transverse isotropy to orthotropy and with 
linear and bimodal sheet models. The results are grouped into 
global and local measures. 

Global measures 

The global measures considered in this study are 1) ejection 
fraction EF, 2) wall thickening WT, and 3) shortening in the long- 
axis direction SL. The variation in EF due to the different fibre 
models is shown in Figure 3(b), with the transversely isotropic 
material law. EF is plotted vs R, and the global fibre models 
SnxRy are grouped according to the exponent n in order to 
highlight the pattern of EF with increasing values of R for each 
value of the exponent ?i = 0.5, 1,2,3. The single value of EF for the 
Karadag fibres is represented as an horizontal line. The predicted 
EF is overall in the range 50%-60%, which can be considered 
physiological [30] , but which is still lower than the measured value 
EF c j m = 77%. Among the global fibre models, for w = 1 ,2,3 
increasing values of the angle R and decreasing values of the 
exponent n raise the value of EF. The reverse holds for n = 0.5, for 
which the increase in the maximum angle R causes a rapid 
decrease of the EF. A possible explanation for this downward 
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trend is provided in the Discussion. The subject-specific Karadag 
model gives one of the highest values of EF together with the 
global models SnlR70 and Sn2R80. 

Changing from the transversely isotropic to the orthotropic 
material law (with the linear sheet model) has the general effect of 
increasing the EF value for all the fibre models, with an average 
increase of 2.7/>p (percentage points), resulting in EF values in the 
range 53%— 63%, i.e. still below the measured value. Figure 4(a) 
shows the comparison of EF between transversely isotropic and 
orthotropic material laws with a linear sheet model for the 
simulations with global fibres with n=\ and Karadag fibres. 
Similar results are obtained for the other values of n. Within the 
orthotropic case, the effect of using a linear or a bimodal sheet 
model has the interesting effect of producing predictions much 
closer to those obtained under transverse isotropy than to 
orthotropy with the linear sheet model, but in general with a 
less consistent pattern across the fibre models (and therefore this is 
not shown). 

For WT, Figure 5(a) represents the wall thickening WT vs R, 
with the transversely isotropic material law. For the global fibre 
models, we observe the same trend as with EF, i.e. for increasing 
values of R and decreasing values of n WT increases. Karadag 
fibres still produce one of the highest values of WT, but some 
global fibres models also give good results. The use of the 
orthotropic material law with the linear sheet model further 
increases the value of WT for every fibre model, as shown in 
Figure 4(b), with an average increase of 43pp. The effect of 
varying the sheet model from linear to bimodal produces again 
values closer to the transversely isotropic case. Overall the range of 
variability for WT is 16%— 37%, with only some fibre models 
giving values close to the physiological range. In any case there is a 
clear discrepancy between the predicted WT values and that 
measured on cine-MRI, WT C j ne = 54%. 

The shortening in the base-to-apex direction (negative values of 
SL) is also substantially affected by the fibre type and material law, 
ranging from — 15% to +10%. Figure 5(b) shows SL vs R, with 
the transversely isotropic material law. Within the global fibre 
models SnxRy, high values of the exponent n combined with low 
values of R produce lengthening (positive values of SL), whereas 
there is increasing shortening for low values of n and increasing 
values of R. In particular, some global fibre models match or 
exceed the measured shortening SL cjne = —10%. Karadag fibres 
produce a shortening value smaller than the measured one. The 
effect of changing form transversely isotropic to orthotropic 
material law (with the linear sheet model) is shown in Figure 4(c). 



The orthotropic material law increases the shortening by 3 . 5pp on 
average. By changing the sheet model to bimodal approximately 
1 pp less shortening was obtained. 

Local measures 

We considered the circumferential strain Ecc, the circumferen- 
tial-longitudinal shear strain Eel and the circumferential-radial 
shear strain Ecr, which were uniformly averaged in the basal 
(upper 35% of the left ventricle), midmyocardial (middle 35% of 
the left ventricle) and apical sections. We show only the average 
strain in the midmyocardial region, because it is a good measure of 
the average behaviour of the left ventricle. In the absence of direct 
strain data, we compared the predicted traces with experimental 
results in the literature. Figure 6(a) shows the Ecc trace for the 
subset of global fibres models with n = 1 and for the Karadag 
model in the case of transverse isotropy. The Karadag model has 
an intermediate behaviour between the fibre models SnlRJO and 
SnlRSO, as we have already found for the global measures. 
Amongst the global fibre models, the magnitude of the strain 
decreases for increasing values of the angle R. All traces start from 
a positive value of Ecc rather than zero, due to the passive 
expansion performed prior to the electromechanical simulation, 
which leaves a residual strain component in Ecc. In order to 
compare predicted Ecc traces with experimental results on Ecc in 
rat, the study by Daire et al. [29], who employed cine-MRI and 
Tagged-MRI, was considered. The measurements of Ecc were 
taken by Daire et al. [29] at fixed percentages of the cardiac cycle, 
therefore it was necessary to convert this temporal subdivision into 
a sequence of time points in the range 0-150/71.? (the duration of 
our simulations). The data (mean and standard deviation) for the 
first half of the cardiac cycle, which corresponds approximately to 
the simulation duration, are shown in red in Figure 6(b). The 
green line is the Ecc trace for the simulation with the Karadag 
fibre model and transverse isotropy, after subtraction of the 
residual strain from the pre-simulation passive expansion. Apart 
from the inevitable inaccuracy of the conversion from cardiac 
cycle percentages to the simulation time interval, there is very 
good agreement between predicted and measured Ecc. 

Figures 7(a) and 8(a) show the Ecr and Eel midventricular traces 
for the subset of global fibres models with n = 1 and for the 
Karadag model in the case of transverse isotropy. With regards to 
Ecr, the Karadag model shows an intermediate behaviour with 
respect to the other fibre models, for which the magnitude of the 
strain decreases for increasing values of the angle R. To the best of 
our knowledge, data on temporal variation of such strains in rat 




Figure 4. Material orthotropy vs transverse isotropy. Effect of introducing sheets by swicthing from transverse isotropy to orthotropy. (a) EF 
(b) WT and (c) SL for global fibre models (with n=\) and Karadag model. 
doi:1 0.1 371 /journal.pone.0092792.g004 
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Figure 5. Global metrics WT and SL. (a) WT and (b) SL for the different fibre models, with transversely isotropic material law. 
doi:1 0.1 371 /journal.pone.0092792.g005 



heart was not available in the literature. Instead, we compared 
predicted shear strains with the measures performed in mouse by 
Zhong et al. [31]. As with the Ecc, it was necessary to convert 
from percentages of cardiac cycle to time points in the range 0- 
\50ms. Figures 7(b) and 8(b) show the comparison between data 
and predicted Ecr and Eel, respectively, for the simulation with 
the Karadag fibre model and transverse isotropy. The agreement 
for Ecr is poorer than for Eel. However, in both cases the trends 
are similar to the experimental ones. 

The effect of sheets was negligible on Ecc, while it was 
substantial on the shear strains. Changing from transverse isotropy 
to orthotropy, and within the orthotropic behaviour, from a linear 
to a bimodal model for the sheets, caused a consistent increase in 
the magnitude of Ecr and Eel, as shown in Figure 9(a) -9(b) for the 
Karadag fibre model in the midmyocardial region. 



Discussion 

In the previous section, results are presented of a systematic 
computational study into the effect of tissue microstructure on 
contraction, using a simplified model of rat left ventricle. It consists 
of a series of electromechanical simulations carried out on a thick- 
walled truncated ellipsoid, with the aim of quantifying the 
combined effect of different models for fibre and sheet directions 
with different tissue passive behaviours during simulated active 
contraction. Sixteen global fibre models were compared against a 
subject-specific model. The impact of sheets was also investigated 
by comparing transversely isotropic and orthotropic material 
properties as well as two models for the sheet orientation, namely a 
linear and a bimodal model. 

In the development of the model we tried to be as species- 
specific as possible in all the components, with the exception of the 




Figure 6. Predicted temporal variation of Ecc. (a) Average midventricular Ecc for global fibre models (with ?1 = 1) and the Karadag model (b) 
Comparison of predicted Ecc for the simulation with the Karadag fibre model and transverse isotropy and measured Ecc in rat [29]. 
doi:1 0.1 371 /journal.pone.0092792.g006 



PLOS ONE | www.plosone.org 



8 



April 2014 | Volume 9 | Issue 4 | e92792 



Effect of Tissue Microstructure on Contraction 




Figure 7. Predicted temporal variation of Ecr. (a) Average Ecr for global fibre models (with n= 1) and the Karadag model.(b) Comparison of 
predicted Ecr for the simulation with the Karadag fibre model and transverse isotropy and measured Ecr in mouse [31]. 
doi:1 0.1 371 /journal.pone.0092792.g007 



reduced range for the cavity pressure, as described in Methods and 
further discussed later on. The geometry has dimensions taken 
from rat cine-MRI scans, the parameters of tissue electrophysiol- 
ogy are based on measurements in rat. Regarding the AP model, 
we preferred the mouse model of Wang and Sobie [21] over other 
available rat models for reasons of computational efficiency. Of 
course, the ventricular activation sequence with simulated apical 
pacing differs from normal in vivo patterns that were present during 
cine-MRI. However, as we compared pre-activation (diastole) with 
peak contraction (systole), we believe that this is unlikely to have 
had a major effect on our findings related to steady-state 
comparisons. As for the mechanics, there are rat-specific 
biophysical contraction models. However, we preferred to use 
the model of Kerckhoffs et al. [24] as it allowed us to adjust its 



parameters in a controlled manner while retaining a good 
description of the contractile behaviour of the myocyte. The 
stiffness parameters of the exponential constitutive law (eqn. (5)) for 
material orthotropy are taken from a study carried out on canine 
hearts by Usyk et al. [4]. In their paper Usyk et al. explore the 
stiffness parameters space in order to study their effect on 
mechanics, providing a range of variation for each of the 
parameters. We chose the orthotropic parameter set that showed 
magnitudes comparable with the values reported for rat in 
transverse isotropy by Omens et al. [32]. This was necessary 
because we did not find published orthotropic stiffness parameters 
values for exponential constitutive laws for rat. We obtained the 
stiffness parameters for transverse isotropy by averaging the 
orthotropic ones in the sheet and normal directions. We ensured 
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that the two parameters sets for orthotropy and transverse isotropy 
had the same magnitudes in order to minimise the effect of 
parameters magnitude and explore only that of material symmetry. 

Various measures were considered, grouped into global and 
local parameters, in order to study the effect of the fibre model, 
material law and sheet model on contraction. For the global 
measures, ejection fraction (EF), wall thickening (WT) and base- 
to-apex shortening {SL) were compared. We chose these measures 
because they are representative of the main features of cardiac 
function in the left ventricle. Furthermore, most electromechanical 
models do not currently predict WT and SL well enough. The 
global measures are consistently affected by the fibre model and 
material law, while they show smaller variation when changing the 
sheet model. The only exception seems to be the downward trend 
of EF in the case of « = 0.5 (Figure 3(b)). We believe this is due to 
an inverse correlation between the magnitude of SL and the rate 
of increase of EF. In other words, the more SL increases in 
magnitude for a given value of n and increasing values of R, the 
smaller the increase in EF. This effect can be seen in Figure 3(b), 
where the steepness of the EF curves decreases consistently going 
from n = 3 to n = 0.5. This observation suggests that there is a 
trade-off between shortening and efficiency of the pumping 
function. Too much shortening is as detrimental to pumping 
function as elongation. Regarding the fibre models, the subject- 
specific fibre model does not produce values of the global measures 
that are notably closer to the measured values than those obtained 
using the best global fibre models. In particular, among the fibre 
models, those with low exponent n and high torsion angle R (such 
as SnlR70 and S11IR8O) produce values of the global metrics that 
are comparable or better than those of the subject specific model. 
These results are encouraging in that they confirm that, in the case 
of the Streeter-like fibre models (see eqn. (1)), the exponent n = \ 
describes fibres in rat left ventricle better than other values, in 
accordance with previous comparisons of rule-based versus DT- 
MRI-derived fibres in rat [33]. In general, the predicted values for 
EF and WT, although being at the physiological range ([30], [13]), 
are lower than the measured ones, while some global fibre models 
are able to match or exceed the measured value for SL. Changing 
from transversely isotropic to orthotropic material properties (with 



the linear sheet model) increases the predicted values of the three 
metrics, but this is not enough to match the measured values for 
EF and WT. Within the orthotropic setting, no significant effect 
was seen on SL, while for EF and WT the bimodal model case 
gave results intermediate between those for transverse isotropy and 
orthotropy with the linear sheet model. 

The Lagrangian strains Ecc, Eel and Ecr are used as local 
measures. The Karadag subject-specific model has strain traces 
close to those of global fibre models with lower exponent n and 
higher angle R. In particular the models SnlR70 and SnlR80 are 
those that in terms of strains give the most similar results to the 
Karadag model. This confirms our findings in the analysis of the 
global metrics. Interestingly, the fibre models that produce higher 
magnitude of strains (that is higher n and lower R) are the same 
that gave poorest results in terms of the global metrics. Conversely, 
this means that fibre models such as Karadag, SnlR70 and 
S11IR8O are able to convert smaller strains into more efficient 
cardiac function. Changing from transversely isotropic to ortho- 
tropic material properties has a negligible effect on Ecc, a small 
effect on Eel, while a more pronounced change was evident in 
Ecr. In particular, introducing the bimodal model increases the 
average magnitude of Ecr in the midmyocardial region by >30% 
with respect to both transverse isotropy and orthotropy with the 
linear sheet model. This is as expected since the effect of 
introducing orthotropy is more likely to have a larger effect on 
Ecr rather than Ecc. In general, the results show strain traces and 
magnitudes close to those reported in the literature for rodents 
[29] [31], especially in the case of Ecc, as shown in Figures 6(b), 
7(b) and 8(b), for Ecc, Ecr and Eel traces, respectively, for the 
simulation with Karadag fibre model and transverse isotropy. 

It is important to discuss the limitations of our model, (a) The 
use of a simplified geometry: the natural geometrical variability of 
the left ventricle is not taken into account in this study, including 
the lack of a neighbouring contracting chamber (right ventricle), or 
of papillary muscles whose orientation is predominantly axial. 
Previous studies suggest this might not bring a significant 
improvement [9], nonetheless, it would be interesting to investi- 
gate if the conclusions about the subject-specific vs global fibre 
models hold also in the case of a subject-specific geometry, (b) The 
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electrical activation pattern: the apical region was chosen for the 
electrical stimulation in order to mimic the apex-to-base direction 
of the electrical activation in the left ventricle. This type of 
activation pattern is consistent, for example, with the type of 
external stimulation employed in the Langendorff-perfused 
experiments. However, this is not what happens in the living 
heart [34]. Evaluating the influence of changes in activation 
pattern on predicted left ventricular mechanics would be an 
interesting application of our model, (c) The ad-hoc procedure to 
estimate the pre-P mesh prior to applying the pressure load: this 
step should be tackled in a more automatic and rigorous approach, 
by solving an inverse problem, (d) The reduced range of pressures 
applied in the model: this is probably the most important 
limitation of the model as pressure is one of the best characterised 
physiological parameters of the heart. If we compare the results of 
the preliminary study we carried out using the model without 
pressure [17] with those presented here, we find similar trends in 
the global measures for the different fibre models. By introducing 
pressure in the model, we find that the magnitude of those 
measures increases significantly on average, meaning that the 
pressure boundary condition is fundamental for the efficiency of 
the simulated pumping function, but does not interfere with the 
effect of the different fibre models on contraction. The use of a 
reduced pressure range is the result of the best compromise 
between realistic contraction patterns and physiological pressures 
in the model. With the current modelling setting, the use of 
realistic (based on laboratory measurements) material properties in 
combination with the chosen contraction model forces the 
pressure to be reduced to avoid non-physiological contraction 
patterns. The alternative would be to enforce physiological 
pressures by modifying the material properties (which would no 
longer be based on experimental measurements) and fully re- 
parametrizing the contraction model (for which experimental data 
would not be available). We chose the first approach for several 
reasons. The first is to reduce the degrees of freedom of the model, 
setting one parameter (the pressure) rather than a set of 
mechanical properties. Another is that of consistency with existing 
literature, which allows us to compare results directiy. Still, as 
mentioned above, this is a major study limitation. As such, it 
identifies a key direction for further study, reinforcing an 
important role of quantitative models as a drivers of scientific 
progress [35], and highlighting that currently available data, 
implemented mechanisms, and attempted explanations are 
insufficient to fully reproduce reality, (e) The lack of the transverse 
angle in the fibre model definitions: this angle regulates the level of 
crossover of fibres which may be important in a correct prediction 
of shear strains [16]. (f) Cross-fibre generation of active tension: at 
the tissue level, the model generates active tension along the fibre 
direction only. Experimental studies have shown that there is a 
certain amount of active tension generated in the cross-fibre plane 
as well [36], and some important electromechanical simulation 
studies already take this into consideration [4], [16]. In particular, 
Bovendeerd et al. [16] have shown that including the cross-fibre 
active tension generation into their model improved the prediction 
of Ecr. (g) The fibre models considered: the choice of rule-based 
and subject-specific fibre models analysed was necessarily reduced 
in order to focus on a specific scientific question. It is important to 
note that we did not consider asymmetrical fibre angle intervals for 
the rule-based model, which have been found experimentally [1] 
and might affect, at least locally, the predicted contraction 
patterns. However, we expect that their effect would remain 
within the ranges predicted by the fibre models already examined. 



Eriksson et al [37], for example, have compared symmetric and 
asymmetric transmural fibre angles for a linear rule-based fibre 
model and found generally consistent predictions in the metrics 
considered. Another important aspect we did not take into account 
is the variability in subject-specific fibre distribution. Karadag et al 
[5] report standard deviations of + 15 degrees for R and + 0.25 
for n. It would be interesting to investigate the inter-subject 
variability of the Karadag fibre model and its effect on contraction 
patterns in the future, (h) The sheet orientation models: we used 
two simplified models of sheet transmural variation, neglecting the 
great regional variability in terms of sheet patterns that has been 
found experimentally, and that is thought to give rise to a 
concertina-like re-arrangement during contraction and relaxation 
[13]. Overall the impact of sheet inclusion in our study was small, 
with the exception of changes in Ecr. It could be argued that this 
could be due to the high stiffness in the sheet direction imposed by 
the constitutive law. However, preliminary tests with different 
magnitudes of the stiffness parameters did not show such a 
substantial variation in the global metrics results to make us believe 
that this problem could be addressed simply by tuning the 
constitutive law parameters. Furthermore, introducing the more 
realistic bimodal sheet model did not affect the global metrics in 
the way we expected. For example, it did not cause a significant 
increase in the predicted wall-thickening, while it has been shown 
that this mechanism is directiy related to sheets dynamics [14]. 

To conclude, we developed and tested an electromechanical 
model and performed a computational study of the impact of 
tissue microstructure on contraction of the left ventricle in rat. We 
showed that the model predicts values of global metrics that are 
within physiological ranges but, except for SL, they fall below the 
values measured on cine-MRI scans. Introducing subject-specific 
fibre models, as opposed to global fibre models, does not provide a 
substantial increase in EF and WT values. The analysis of strains 
confirmed that the subject-specific fibre model gave predictions 
similar to linear global fibre models with higher rotation angle. 
Changing from transverse isotropy to orthotropy in the passive 
behaviour, therefore introducing the effect of tissue anisotropy due 
to the sheets, improves the predictions of EF and WT, but they 
still remain lower than the measured values. This means that 
including the sheet orientation in the model is an important factor, 
and a more detailed representation of angular re-orientation of 
sheet material may be required. Therefore, we believe that, within 
the current modelling framework, the subject-specific approach in 
fibre definition by itself is not able to compensate for the model 
limitations that are believed to come from the underestimation of 
sheet dynamics, and that further investigation of the dynamics of 
sheet orientation should be carried out and novel approaches 
developed in order to fully represent its effect on cardiac 
mechanics. 
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